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ABSTRACT 

We consider the single-spin-flip dynamics of the random-field Ising model on a 
Bethe lattice at zero temperature in the presence of a uniform external field. We 
determine the average magnetization as the external field is varied from — cxd to 
+00 by setting up the self-consistent field equations, which we show are exact in 
this case. The qualitative behavior of magnetization as a function of the external 
field unexpectedly depends on the coordination number z of the Bethe lattice. 
For z = 3, with a gaussian distribution of the quenched random fields, we find 
no jump in magnetization for any non-zero strength of disorder. For z > 4, 
for weak disorder the magnetization shows a jump discontinuity as a function of 
the external uniform field, which disappears for a larger variance of the quenched 
field. We determine exactly the critical point separating smooth hysteresis curves 
from those with a jump. We have checked our results by Monte Carlo simulations 
of the model on 3- and 4- coordinated random graphs, which for large system sizes 
give the same results as on the Bethe lattice, but avoid surface effects altogether. 
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I. Introduction 



Recently, a simple model has been introduced [1] for hysteresis in magnets, 
which incorporates interesting effects like the return-point memory and Barkhausen 
noise [2]. In this model, Ising spins with a quenched random field at each site 
evolve by a zero-temperature single-spin-flip dynamics. The authors argued that 
in this model, if the external field is increased slowly, the steady-state magneti- 
zation as a function of the field has a jump discontinuity at some critical value of 
the field for a small disorder, but is a continuous function with no jump discon- 
tinuity for large disorder. This picture was supported by numerical simulations 
of the model on hypercubic lattices in two and higher dimensions. Subsequent 
work [3] studied in detail the transition from jump to no-jump in magnetization 
at a critical value of the gaussian disorder, and observed scaling behavior in the 
neighbourhood of this critical disorder. However, the exact solution of this model 
in one dimension does not show a jump discontinuity for gaussian disorder neither 
for the ferromagnetic nor for antiferromagnetic exchange couplings [4]. 

In this paper, we extend the treatment of [4] to study this hysteresis model on 
a Bethe lattice. The Bethe lattice of coordination z is the formal infinite-size limit 
of a branching tree (the Cayley tree) where each spin has z nearest neighbors, and 
the statistical averages are calculated away from the 'surface'of the lattice. For 
^ = 3, the lowest non-trivial coordination, we find no jump in magnetization for 
any nonzero disorder, if the quenched random fields have a gaussian distribution 
— just as in the one-dimensional case. For coordination 2; > 4, we find there is 
a non-zero critical disorder where the macroscopic jump in magnetization in the 
hysteresis loop first disappears. 

This is very surprising, as in all the models studied on the Bethe lattice, the 
qualitative behavior of the solution has been found to be independent of the 
coordination number ( so long as it is greater than 2). In particular, a Bethe 
lattice with finite coordination number z has the same critical behavior as the 
mean-field theory, which corresponds to the limit of large coordination number z 
and coupling constant scaling as Xjz. The reason why this unusual dependence 
on z shows up in this problem is not yet understood. 

Our treatment is based on setting up self-consistent equations for some near- 
est neighbour correlation function in the problem [5]. We can show that these 
self-consistent equations are exact in this case, though we are not aware of a 
rigorous proof that this happens in general for a Bethe lattice in the presence 
of quenched disorder. In fact, the presence of the disorder usually renders the 
problem analytically intractable. For example, for the Ising spin-glass problem 
on a Bethe lattice with random ± J bonds, it has not been possible to determine 
exactly even the zero temperature quantities like the ground state energy or the 
ground state entropy [6,7]. However, the qualitative behavior of the system near 
the thermal critical point seems fairly well understood [8]. 
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The plan of this paper is as follows: In Section II, we define the model pre- 
cisely. In Section III, we set up recursion relations on a Cayley tree for conditional 
probabilities that the spin at a given site at height r from the boundary is down, 
given that the spin on its parent, "upward" neighbor on the tree is down. We 
are interested in the intensive quantities, such as magnetization or energy density 
on the tree far away from the boundary. These turn out to be independent of 
details of the boundary conditions, and we take this as the definition of Bethe 
approximation in our case. We obtain an explicit expression for magnetization 
as a function of external field for arbitrary distribution of the quenched random 
fields. In Section IV, we describe a method to simulate spin systems on the Bethe 
lattice that is computationally efficient, and does not suffer from surface effects. 
We use this method to check the validity of our self-consistent equations for the 
case of gaussian and rectangular distributions of quenched random fields. The 
agreement is found to be excellent. Section V contains some concluding remarks. 



II. The Model 

We consider a lattice of N sites. Each site is labeled by an integer i — 
1 to N, and carries an Ising spin Si {Si — ±1) which interacts with a finite 
number z of neighbouring spins with a ferromagnetic interaction J. There is 
a uniform magnetic field h which is applied externally. In addition, at each 
site i, there is a local quenched random field hi. The fields {hi} are assumed 
to be independent identically distributed random variables with a continuous 
probability distribution p{h). The system is described by the Hamiltonian 

h^-jY. s,Sj -Y^hiSi-hY, Si. (1) 

(ij) i i 

The zero-temperature single-spin-fiip Metropolis-Glauber dynamics [9] is speci- 
fied by the transition rates 

Rate [Si -Si] = T, ff l^E < 0; 

= 0, otherwise; (2) 

where AE is the change of energy of the system as a result of the spin-fiip. We 
shall be interested in long time scales ^ F"^. In this limit, the dynamical rule 
simplifies to the following: Choose a spin at random, and fiip it only if this process 
would lower the energy. Repeat till a stable configuration is obtained. 

The problem of hysteresis which we address here is as follows: Start with a 
sufficiently large negative applied field h, so that in the stable configuration all 
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spins are down {Si — —1, for all i) and increase the field slowly. At some value 
of h, the local field ii at some site i, defined by 



J'^Sj + hi + h 



(3) 



j 



will become positive, and this spin would flip up. [The summation in (3) is over 
all the neighbors j of i]. This changes the effective field at the neighbors, and 
some of them may flip up, and so on, causing an avalanche of flipped spins. We 
determine the total magnetization when the avalanche has stopped. Then we 
raise the applied field a bit more, and determine the magnetization in the stable 
state again. The process is continued until all the spins flip up. This generates 
the lower half of the hysteresis loop (plot of magnetization m{h) versus h) in the 
situation where the applied field is varied very slowly, or equivalently, when the 
spins relax infinitely fast. The upper half of the hysteresis loop is obtained 

when the field h is decreased from +00 to —00. This is related to the lower half 
of the loop mi{h) by symmetry 



This corresponds to the zero frequency limit of a driving field oscillating si- 
nusoidally in time with frequency uj. Note that the limit of a; —> is taken after 
the limit temperature T — > 0. If the limits are taken in the reverse order, the 
area of the hysteresis loop goes to zero as u goes to zero for all non-zero T. 

An important feature of the above dynamics for ferromagnetic couplings is 
that if we start with any stable configuration, and then increase the external field 
and allow the system to relax, then in the relaxation process no spin fiips more 
than once. Furthermore, the final stable configuration is the same whatever the 
order in which unstable spins are fiipped. This property is called the 'no passing 
property' [10], and greatly simplifies the analysis. 

III. Recursion Relations on the Cayley Tree 

The standard approach for solving statistical mechanics problems on the 
Bethe lattice is to consider the problem on a Cayley tree, and consider behavior 
deep inside the tree, i.e. far from the boundaries of the tree [11,12]. If suit- 
able care is taken to remove the effects of the boundary, all correlation functions 
deep inside the Cayley tree (say for the Ising model with an external field) are 
found to be the same as in the Bethe approximation. Thus, we may say that the 
Bethe lattice is the deep interior part of the Cayley tree. Here we shall use this 
approach. In Section IV, a different approach is presented. 

Consider a Cayley tree of height n. Each site of the tree has coordination 
number z, except the boundary sites which have coordination number 1. The 



m„(/i) 



mi{-h). 



(4) 
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level n consists of only one site O, called the central site. For r > 1 the level 
(n — r) has exactly z ■ {z — 1Y~^ sites [Fig. 1]. 

Wc start with the external field h large and negative, so that ground state 
of the system is with all spins down. Now, increase the external field to a finite 
value /i, and flip up any spin for which the net local field is positive. As the 
same final stable configuration is attained, whatever the order in which spins are 
relaxed, we may start by first relaxing spins of level 1. Then we relax spins of 
level 2, then of level 3, and so on. If a spin at level r is fiipped up, we check all 
its descendents again for possible upward fiips. 

Let Pr be the conditional probability that a randomly chosen spin at level r is 
upturned in this scheme, given that its parent spin at level (r + 1) is kept down, 
and the spin and all its descendent spins are relaxed as far as possible. Let 
be the spin at level r. We relax all the descendent spins of S*,. first, keeping S.^ 
down. In this process, each of the z — 1 direct descendents of at level (r — 1) 
is independently flipped up with probability P^-i- Hence the probabilities that 
2; — 1,2; — 2,.. .,0 of the children of Sj. are flipped up in this relaxation process are 
Pr-i, - l)P^^rf (1 - Pr-i), . . . , (1 - Pr-iY'^ respectively. Consider the case 
where s of the children are up: since the parent neighbor remains down for this 
part of the calculation, the net number of down neighbors is z — 2s, and hence, 
the spin Sr will flip up if the local field at this site exceeds {z — 2s)J — h. Let 
Ps{h) denote the probabihty that the local field at a randomly chosen site is large 
enough so that the spin will flip up if s of its children are up, and the uniform 
held is h. Clearly 

Ps{h) — Prob that local held > —h + zJ — 2s J 

/oo 
p{hi)dhi. (5) 
-h+{z-2s)J 

Then it is easily seen, e.g. ior z — 3 that 

= 3) : Pr = P^.Mh) + 2P.-i(l - Pr-i)pi{h) + (1 - Pr-ifpo{h). (6) 

Given a value of h, we determine the quantities Ps{h). Then using Eq. (6), and 
the initial condition Pi = pi{h), we can recursively determine Pr for all r > 2. 
For large r ^ n, Pr tend to a flxed point P* given by the self-consistent equation 

= E (r') ^"(1 - py-'-^'Prih) (7) 

r=0 

This is a polynomial equation in P*, which can be solved in terms of {ps{h)}. 
Finally, for the central site O at level n, there are z children, and a similar 
argument gives 

Prob(5o = +1) = E ir) P'-'il - Py-'Pr{h) (8) 

r=0 
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Substituting the value of P*, from Eq. (7), we determine the probabihty that 
this spin Sq is up, and hence the average magnetization at this site. 

The arguments above do not require that all the z descendent subtrees of O 
be of equal height. So long as O is sufficiently far from the boundary, we get the 
same conditional probability P*, and hence the same value of magnetization. This 
proves that all sites 'deep inside' the tree have the same average magnetization. 

IV. Simulations 

The derivation of our self-consistent equations assumes the existence of a 
unique thermodynamic state deep within the Cayley tree which is independent 
of boundary conditions. While this is quite plausible, uniqueness of the Gibbs 
state has been proved so far for the RFIM only in one dimension, and only for 
a bivariate distribution of the quenched field, and nonzero temperatures [13]. It 
seems desirable to have a direct check of these equations by numerical simulations 
which do not involve making any assumptions about the thermodynamic state. 

While the procedure of the previous section treating the Bethe lattice as 
sites deep inside the Cayley tree is well known and conceptually simple, it is 
not suited for numerical simulations. Most of the sites of the Cayley tree are 
within a short distance from the surface, and cannot be used for averaging. Since 
the 'bulk' forms a negligible fraction of all possible sites, special care has to 
be taken to subtract the surface contribution. For our simulations, we used a 
different technique that is computationally efficient and gets rid of surface effects 
altogether. This technique has been used earlier to study spin systems on random 
graphs by Monte Carlo simulations [14]. 

Our simulation algorithm involves construction of a random graph having N 
sites such that each site has exactly z neighbors. The precise algorithm we used 
was as follows: Label the N sites by integers from 1 to N. We shall assume is 
even in the following. Connect site i to site (i + 1) for all i. Site N is connected 
to site 1. This gives us a ring of N sites. Now construct {z — 2) independent 
random pairing of N sites into N/2 pairs, and add a bond for each of the paired 
sites. Thus, we get a graph in which each site has coordination number z [Fig. 

2]. 

In this construction, all sites are on same footing, and there is no 'surface'. 
Unlike Cayley tree, this graph has loops. However it is easy to see that there are 
typically very few small loops. For example, for z = 3 the probability that sites 
i, {i + 1) and {i + 2) form a loop of length 3 is the probability that site i is paired 
with {i + 2), and equals 1/{N — 1). Thus the expected number of loops of size 
3 in a 2; = 3 graph of N sites tends to 1 for large N. Similarly, it can be shown 
that the expected number of loops of length 4 is 2 for large N. In general, the 
average number of loops of length £ increases as with A = 2; — 1 for the random 
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graph with coordination number z, and is a neghgible fraction of all sites belong 
to any loop of length < £ lor £ <^ log A^/ log A [15]. 

If the smallest loop going through a given site is of length < {2d +1), then it 
follows that up to a distance d from that site, the lattice looks like a Bethe lattice. 
Hence our random lattice would look like a Bethe lattice for 2; = 3 at almost all 
sites for a distance <, log2 N . This, in turn, can be shown to imply that in the 
thermodynamic limit N — > cxo, the free energy per site on our lattice for classical 
statistical mechanical models with short range interactions (say nearest neighbor 
only) are the same as in the Bethe-Peierls approximation. 

In our simulations, we used N — 10^. We used simple scanning to decide 
which spins to be flipped at the next time step. The dotted hues in Figs. 3 
and 4 show the results of a simulation for ;i; = 3 for quenched gaussian random 
fields with mean and variance a = 1 and (7 = 3 respectively. The lower and 
upper halves of the hysteresis loop were obtained separately in the simulation. 
Also shown in the figures are the results of solution of Eqs. (7-8). The statistical 
errors of the simulation are quite small. Different runs, with different realizations 
of quenched fields give results which are indistinguishable at the scale of the 
graph. The agreement with the theoretical calculation is excellent. For much 
smaller values of disorder cr <, .1, the hysteresis loops are very approximately 
rectangular. In this case, the value of coercive field is governed by the largest 
realized value of quenched local field, which shows noticeable sample to sample 
fluctuations. As noted above, for z — the hysteresis loop is smooth for all values 
of the disorder greater than zero: the quadratic equation (6) has only one stable 
solution. 

For 2; = 4, we do find a transition. At small disorder, the hysteresis loop 
has a jump: one large event fiips a large, finite fraction of the spins in the 
thermodynamic limit. Figure 5 shows the analytical and simulation results for 
a = 1.75; there is a large jump in the magnetization at /i = 1.0037. The critical 
value of the disorder is very slightly larger than this (cr^^^^-* = 1.78126) and above 
the critical disorder the hysteresis loop is smooth. The critical field he — 1 ai 
the critical disorder for 2; = 4 and we conjecture for larger z as well. This simple 
result follows from the observation that at h = 1, P = 1/2 is always a fixed 
point for Eq. (7) for all z: (Tc{z) may be determined by making this fixed point a 
double root. This makes the transition a traditional saddle-node transition (the 
lower branch merges with an "unstable" branch of the self-consistent equation). 
The critical exponents are thus the same as that for the infinite-range mean- 
field model (which also undergoes a saddle-node transition in its self-consistent 
equation). We have checked that this same pattern also occurs for z = 5 (where 
it gives o"c = 2.58201), and conjecture that it gives the correct critical point for 
all 2 > 3; in 2 = 3 the coalescence between the stable and unstable branches of 
the M{h) curves never occurs. 

V. Discussion 
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It is natural to compare the zero temperature hysteresis on the Bethe lattice 
with the corresponding infinite-range mean field result obtained in the limit of 
large coordination number when the ferromagnetic coupling is taken to be J /N , 
the same for all pairs of sites. In this case the mean-field solution is given by [1] 



For a < Gc = the above equation has two solutions m\{h) and m*^{h) 

which are related to each other by the symmetry m}{—h) = —m*^{h). These 
correspond to the two halves of the hysteresis loop for increasing and decreasing 
field respectively. 

For a > (Tc, Eq. (9) has a single valued real solution m*(h) which is an odd 
function of h. Thus there is no hysteresis for a > a^- The remanence goes to zero 
continuously as a tends to a"c from below. For a < (7^ there is a discontinuity 
in magnetization at a critical field he- The value of he, and the magnitude of 
jump in magnetization both tend to zero continuously as o" — > cTc- This lack of 
hysteresis for cr > dc is an artifact of the hard-spin infinite-range model: our 
Bethe lattice has hysteresis at all values of a (as does the infinite-range model 
with continuous spins in a double- well potential [1]). 

For z > 3 the behavior as one approaches Uc from below is similar to the 
infinite-range model, and the corresponding critical exponents (3 and 5 will be 
the same [1]. 

So far, we have discussed the case when the quenched random fields have 
an unbounded distribution. For bounded distributions of the quenched random 
field, one can get jumps in magnetization even for z — 2>. Consider, for example, 
the case when {hi} have a uniform rectangular distribution between —hmax and 
+/im,aa;[16]. If wc Start with all spins down, and increase field slowly, clearly 
nothing happens for h < 3J — h^ax- If h exceeds this value, then the spin with 
the largest value of hi will flip up. If hmax < J: this will make the net local field at 
the neighbors positive. These spins will flip up, which in turn fiips their neighbors, 
and so on. Thus, for h^ax < J, the magnetization m[h) jumps discontinuously 
from —1 to +1, as h cross 3J — hmax [17]. If J < hmax < 2J, one can show 
that same thing occurs as the system is, on the average, unstable for creating 
such a 'nucleus' of up spins. However, for hmax > 2J, this particular instability 
is absent, and the magnetization is a continuous function of h. Note that the 
magnetization jump goes discontinuously to zero. If the distribution of quenched 
fields p{hi) has delta functions, in addition to a continuous part, clearly this will 
lead to discontinuities in the m{h) curve. Any other singularities oip{hi), say at 
hi = a lead to singularities in m{h) for /i = a ± 3 J, a ± J. 

An interesting open question, which we have not been able to answer so far 
is to characterize all possible 'metastable' states on the Bethe lattice. Are all 
of these obtainable as solutions of self-consistent equations of the type discussed 



m = erf 



Jm + h 



(9) 



V2^ 
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above? For example, can one calculate the magnetization when the external field 
is first increased monotonically from —oo to a value Hi, and then reduced to 
a value H2 < Hil Further study of such questions would perhaps help in our 
understanding of the more general question of hysteretic dynamics of systems 
with many metastable states. 
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Captions to figures 



Fig. 1: A Cayley tree of coordination number 3 and height 3. 

Fig. 2: An example of a random graph with coordination number 3. Dotted hues 
indicate the random pairs. 

Fig. 3: Hysteresis loop on the Bethe lattice of coordination number 3. Case shown 
is for standard deviation of quenched random field a — J. The result 
of simulation for = 10® spins (points) is in good agreement with our 
theoretical result (continuous curve). 

Fig. 4: Hysteresis loop on the Bethe lattice for 2; = 3 and a — 3J. The result 
of simulation for A^ = 10® spins (points) is in good agreement with our 
theoretical result (continuous curve). 

Fig. 5: Magnetization curves for the Bethe lattice of coordination number 4 in 
increasing field. 
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